Active Contours using Level Sets

This tour explores image segementation using level set methods.

Installing toolboxes and setting up the path.

You need to download the following files: signal toolbox, general toolbox and graph toolbox.

You need to unzip these toolboxes in your working directory, so that you have toolbox_signal, toolbox_general and toolbox_graph in your directory.

For Scilab user: you must replace the Matlab comment '%' by its Scilab counterpart '//'.

Recommandation: You should create a text file named for instance numericaltour.sce (in Scilab) or numericaltour.m (in Matlab) to write all the Scilab/Matlab command you want to execute. Then, simply run exec('numericaltour.sce'); (in Scilab) or numericaltour; (in Matlab) to run the commands.

Execute this line only if you are using Matlab.

getd = @(p)path(p,path); % scilab users must *not* execute this

Then you can add the toolboxes to the path.

getd('toolbox_signal/');
getd('toolbox_general/');
getd('toolbox_graph/');

Managing level set functions

In order to perform curve evolution, we will deal with a distance function stored in a 2D image D. The curve is embedded in the level set D=0. This curve will be evolved by modifying the image D. A curve evolution ODE can be replaced by an PDE on D. This allows to deal with topological changes when a curve splits or two curves merge.

Size of the image.

n = 200;
[Y,X] = meshgrid(1:n,1:n);

First load a cirle shape.

% radius
r = n/3;
% center
c = [r r] + 10;
% shape
Dcirc = sqrt( (X-c(1)).^2 + (Y-c(2)).^2 ) - r;

Exercice 1: (the solution is exo1.m) Load a square shape Dsq at a different position for the center, e.g. n - 10 - [r r].

exo1;

Display the curves.

clf;
subplot(1,2,1);
plot_levelset(Dcirc);
subplot(1,2,2);
plot_levelset(Dsq);

Exercice 2: (the solution is exo2.m) Compute the intersection and the union of the two shapes. Store the union in D0, that we will use in the remaining part of the tour.

exo2;

Mean Curvature Motion.

The mean curvature motion tries to shorten the length of the level sets of a function. In particular it shrinks the object represented here by the level set 0.

The mean curvature motion of the level sets of some image D is driven be the following equation.

 |dD/dt = norm(grad(D)) * div(grad(D)/norm(grad(D)) ) |

Parameters for the descent of the energy.

% maximum time of evolution
Tmax = 200;
% time step (should be small)
tau = .5;
% number of iterations
niter = round(Tmax/tau);
% use centered differences
options.order = 2;
% initial shape at t=0
D = D0;

Compute the right hand side of the evolution equation.

% the gradient
g0 = grad(D,options);
% norm of the gradient
d = max(eps, sqrt(sum(g0.^2,3)) );
% normalized gradient
g = g0 ./ repmat( d, [1 1 2] );
% the curvature
K = d .* div( g,options );

Perform the descent (this correspond to an explicit discretization in time of the PDE).

D = D + tau*K;

Exercice 3: (the solution is exo3.m) Implement the mean curvature motion.

exo3;

Levelset Re-distancing

During PDE resolution, a level set function D might become ill-conditionned, so that the zero crossing is not sharp enough. The quality of the level set function is restored by computing the signed distance function to the curve D==0.

We simulate a decrease of quality by raising to a power 3 the distance function.

D = D0.^3;

We compute the distance function.

D1 = perform_redistancing(D0);

Display the level sets.

clf;
subplot(1,2,1);
plot_levelset(D);
title('Before redistancing');
subplot(1,2,2);
plot_levelset(D1);
title('After redistancing');

Edge-based Segmentation with Geodesic Active Contour

Geodesic active contours are basically snakes and level set in the same framework.

Note: these active contours should not be confounded with the geodesic shortest paths, that are geodesics between two points. Here the active contour is a close curve progressively decreasing a weighted geodesic length.

First we load an image to segment.

n = 200;
name = 'cortex';
M = rescale( sum( load_image(name, n), 3) );

Given a background image M to segment, one needs to compute an edge-stopping function E. It should be small in area of high gradient, and high in area of large gradient.

Exercice 4: (the solution is exo4.m) Compute an edge-stopping function E by smoothing the inversed magnitude of the gradient 1./(epsilon+norm(grad(M)). Rescale E so that it ranges [.1,1].

exo4;

The geodesic active contour minimizes the length of curve.

min_D \int_{D==0} E

The geodesic active contour evolution is a descent of this energy, and corresponds to a mean curvature motion modulated by the edge-stopping function.

Exercice 5: (the solution is exo5.m) Compute an initial shape, for instance a square centered at [n n]/2.

exo5;

We set parameters for the descent.

% final time
Tmax = 1500;
% time step
tau = .4;
% number of steps
niter = round(Tmax/tau);
% initial distance function
D = D0;

The gradient of the geodesic active contour energy reads:

|dD/dt = E * norm(grad(D)) * div( grad(D) / norm(grad(D)) ) + grad(E),grad(D) |

Exercice 6: (the solution is exo6.m) Compute this gradient G (right hand side of the PDE) using the current value of the distance function D.

exo6;

Do the descent.

D = D + tau*G;

Once in a while (e.g. every 30 iterations), perform re-distancing.

D = perform_redistancing(D);

Exercice 7: (the solution is exo7.m) Implement the geodesic active contours gradient descent. Do not forget to do the re-distancing.

exo7;

Region-based Segmentation with Chan-Vese

Chan-Vese active contours are basically Mumord-Shah + level sets.

The geodesic active contour uses an edge-based energy. It has lots of local minima and is very sensitive to initialization. In order to circumvent these drawbacks, one can use a region based energy like the Mumford-Shah functional. Re-casted into the level set framework, it reads:

 |min_D Area(D>0) + lambda*\int_{D>0} |M-c1| + lambda*\int_{D<0} |M-c2|^2|

Where c1 and c2 are the expexted values

Exercice 8: (the solution is exo8.m) Compute an initial shape, for instance many small circles.

exo8;

We set parameters for the energy.

lambda = 0.8;
c1 = 0.7; % black
c2 = 0; % gray

We set parameters for the descent.

% final time
Tmax = 100;
% time step
tau = .4;
% number of steps
niter = round(Tmax/tau);
% initial distance function
D = D0;

The corresponding gradient descent is the Chan-Vese active contour method:

dD/dt = norm(grad(D))*( div(grad(D)/norm(grad(D))) - lambda*(M-c1)^2 + lambda*(M-c2)^2 )

Exercice 9: (the solution is exo9.m) Compute this gradient G using the current value of the distance function D.

exo9;